/*
 * solverThread_cpu.c
 *
 *  Created on: 2011-3-17
 *      Author: zhang
 */
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <kpm_plans.h>
#include <timer.h>
#include <mtrnd.h>
#include <ell.h>
#include <randomseed.hpp>

#include <math.h>

template<typename T>
inline T createGussianRand(MT::MersenneTwist &mt)
{
  return sqrt(-2 * log(mt.genrand_real3())) * cos(2 * M_PI * mt.genrand_real3());
}

template<typename T>
T* solveDOSPlan_cpu(DOSPlan<T>* plan)
{
	MT::MersenneTwist mt;
	mt.init_genrand(RandomSeed::generate());

	ELL<T>* H = plan->H;
	int h_size = H->rowsNum;
	int rs = plan->num_r;
	int M = plan->m;
	memset(plan->U, 0 ,sizeof(T) * M);

	T* VR = (T*)malloc(sizeof(T) * h_size);
	T* TempR1 = (T*)malloc(sizeof(T) * h_size);
	T* TempR2 = (T*)malloc(sizeof(T) * h_size);

	if(VR == NULL || TempR1 == NULL || TempR2 == NULL)
	{
        printf("failed to malloc memory for vectors.");
        exit(0);
	}

	for(int s=0; s<rs; s++)
	{
	    //generate R0 and perform vector-vector multiplication.
	    for(int j=0;j<h_size;j++)
	    {
	      TempR1[j] = VR[j] = createGussianRand<T>(mt);
	      plan->U[0] += VR[j] * VR[j];
	    }

		//generate R1 and perform vector-vector multiplication.
		for(int i=0;i<H->rowsNum;i++){
		  T sum = 0.0;
		  int k = i * H->colsNum;

		  for(int j=k; j<k+H->colsNum; j++){
		  	 sum += H->elements[j]*TempR1[H->cols[j]];
		  }


/*
		  sum += H->elements[k]*TempR1[H->cols[k]];
		  sum += H->elements[k+1]*TempR1[H->cols[k+1]];
		  sum += H->elements[k+2]*TempR1[H->cols[k+2]];
		  sum += H->elements[k+3]*TempR1[H->cols[k+3]];
		  sum += H->elements[k+4]*TempR1[H->cols[k+4]];
		  sum += H->elements[k+5]*TempR1[H->cols[k+5]];
		  sum += H->elements[k+6]*TempR1[H->cols[k+6]];
*/
		  TempR2[i] = sum;
		  plan->U[1] += VR[i] * TempR2[i];
		}

		for(int n=2;n<M;n++)
		{

			for(int i=0;i<H->rowsNum;i++){
			  int k= i * H->colsNum;
			  T sum = 0.0;
			  for(int j=k; j<k+H->colsNum; j++){
				 sum += H->elements[j]*TempR2[H->cols[j]];
			  }
/*
			  sum += H->elements[k]*TempR2[H->cols[k]];
			  sum += H->elements[k+1]*TempR2[H->cols[k+1]];
			  sum += H->elements[k+2]*TempR2[H->cols[k+2]];
			  sum += H->elements[k+3]*TempR2[H->cols[k+3]];
			  sum += H->elements[k+4]*TempR2[H->cols[k+4]];
			  sum += H->elements[k+5]*TempR2[H->cols[k+5]];
			  sum += H->elements[k+6]*TempR2[H->cols[k+6]];
*/
			  TempR1[i] = sum * 2 - TempR1[i];
			  plan->U[n] += VR[i] * TempR1[i];
			}

		   T* tp = TempR1;
		   TempR1 = TempR2;
		   TempR2 = tp;
		}
	}

	for(int i=0;i<M;i++){
		plan->U[i] /= (T)(rs * h_size);
		printf("%.3f ", plan->U[i]);
	}

	free(VR);
	free(TempR1);
	free(TempR2);
	return 0;
}

template double* solveDOSPlan_cpu<double>(DOSPlan<double>* plan);
template float* solveDOSPlan_cpu<float>(DOSPlan<float>* plan);
